#include <stdio.h>
#include <string.h>
#include <stdint.h>
#include <utils.h>
#include <sm2_bn.h>


const sm2_bn_type SM2_BN_ONE = {1, 0, 0, 0};

const uint64_t *sm2_bn_one(void) {
    return &SM2_BN_ONE[0];
}

void sm2_bn_set_zero(sm2_bn_type r) {
    r[0] = 0, r[1] = 0, r[2] = 0, r[3] = 0;
}

int sm2_bn_rand_range(sm2_bn_type r, const sm2_bn_type range) {
    unsigned int tries = 100;
    do {
        if (!tries) {
            return 0;
        }
        tries--;
        if (gen_random_bytes((uint8_t *) r, 32) != 1) {
            print_log();
            return -1;
        }
    } while (sm2_bn_compare(r, range) >= 0);
    return 1;
}

void sm2_bn_from_bytes(sm2_bn_type r, const uint8_t in[32]) {
    r[0] = GET_U64(in + 24), r[1] = GET_U64(in + 16), r[2] = GET_U64(in + 8), r[3] = GET_U64(in);
}

void sm2_bn_to_bytes(const sm2_bn_type a, uint8_t out[32]) {
    PUT_U64(out, a[3]), PUT_U64(out + 8, a[2]), PUT_U64(out + 16, a[1]), PUT_U64(out + 24, a[0]);
}

void sm2_bn_copy(sm2_bn_type d, const sm2_bn_type s) {
    d[0] = s[0], d[1] = s[1], d[2] = s[2], d[3] = s[3];
}

void sm2_bn_copy_conditional(sm2_bn_type dst, const sm2_bn_type src, uint64_t move) {
    const uint64_t mask_one = 0 - move;
    const uint64_t mask_two = ~mask_one;
    dst[0] = (src[0] & mask_one) ^ (dst[0] & mask_two);
    dst[1] = (src[1] & mask_one) ^ (dst[1] & mask_two);
    dst[2] = (src[2] & mask_one) ^ (dst[2] & mask_two);
    dst[3] = (src[3] & mask_one) ^ (dst[3] & mask_two);
}

static uint64_t uint64_t_is_zero(uint64_t in) {
    in |= 0 - in;
    in = ~in;
    in >>= 63;
    return in;
}

uint64_t sm2_bn_is_equal(const sm2_bn_type a, const sm2_bn_type b) {
    const uint64_t res = (a[0] ^ b[0]) | (a[1] ^ b[1]) | (a[2] ^ b[2]) | (a[3] ^ b[3]);
    return uint64_t_is_zero(res);
}

int sm2_bn_compare(const sm2_bn_type a, const sm2_bn_type b) {
    if (a[3] < b[3]) return -1;
    if (a[3] > b[3]) return 1;
    if (a[2] < b[2]) return -1;
    if (a[2] > b[2]) return 1;
    if (a[1] < b[1]) return -1;
    if (a[1] > b[1]) return 1;
    if (a[0] < b[0]) return -1;
    if (a[0] > b[0]) return 1;
    return 0;
}

uint64_t sm2_bn_is_zero(const sm2_bn_type a) {
    return uint64_t_is_zero(a[0]) & uint64_t_is_zero(a[1]) & uint64_t_is_zero(a[2]) & uint64_t_is_zero(a[3]);
}

uint64_t sm2_bn_plus(sm2_bn_type result, const sm2_bn_type a, const sm2_bn_type b) {
    uint64_t carry = 0, temp = a[0] + b[0];
    carry = temp < a[0];
    result[0] = temp;

    temp = a[1] + carry;
    carry = temp < a[1];
    result[1] = temp + b[1];
    carry += result[1] < temp;

    temp = a[2] + carry;
    carry = temp < a[2];
    result[2] = temp + b[2];
    carry += result[2] < temp;

    temp = a[3] + carry;
    carry = temp < a[3];
    result[3] = temp + b[3];
    carry += result[3] < temp;
    return carry;
}

uint64_t sm2_bn_subtraction(sm2_bn_type result, const sm2_bn_type a, const sm2_bn_type b) {
    uint64_t carry = 0, temp = a[0] - b[0];
    carry = temp > a[0];
    result[0] = temp;

    temp = a[1] - carry;
    carry = temp > a[1];
    result[1] = temp - b[1];
    carry += result[1] > temp;

    temp = a[2] - carry;
    carry = temp > a[2];
    result[2] = temp - b[2];
    carry += result[2] > temp;

    temp = a[3] - carry;
    carry = temp > a[3];
    result[3] = temp - b[3];
    carry += result[3] > temp;

    return carry;
}

void sm2_bn_multiplication(sm2_bbn_type result, const sm2_bn_type a, const sm2_bn_type b) {
    uint64_t a_[8], b_[8], s[16] = {0};

    for (int i = 0; i < 4; i++) {
        a_[2 * i] = a[i] & 0xffffffff;
        b_[2 * i] = b[i] & 0xffffffff;
        a_[2 * i + 1] = a[i] >> 32;
        b_[2 * i + 1] = b[i] >> 32;
    }

    for (int i = 0; i < 8; i++) {
        uint64_t u = 0;
        for (int j = 0; j < 8; j++) {
            u = s[i + j] + a_[i] * b_[j] + u;
            s[i + j] = u & 0xffffffff;
            u >>= 32;
        }
        s[i + 8] = u;
    }

    for (int i = 0; i < 8; i++) {
        result[i] = (s[2 * i + 1] << 32) | s[2 * i];
    }
}

static uint64_t sm2_bbn_plus(sm2_bbn_type r, const sm2_bbn_type a, const sm2_bbn_type b) {
    uint64_t carry = 0, temp = a[0] + b[0];
    carry = temp < a[0];
    r[0] = temp;

    temp = a[1] + carry;
    carry = temp < a[1];
    r[1] = temp + b[1];
    carry += r[1] < temp;

    temp = a[2] + carry;
    carry = temp < a[2];
    r[2] = temp + b[2];
    carry += r[2] < temp;

    temp = a[3] + carry;
    carry = temp < a[3];
    r[3] = temp + b[3];
    carry += r[3] < temp;

    temp = a[4] + carry;
    carry = temp < a[4];
    r[4] = temp + b[4];
    carry += r[4] < temp;

    temp = a[5] + carry;
    carry = temp < a[5];
    r[5] = temp + b[5];
    carry += r[5] < temp;

    temp = a[6] + carry;
    carry = temp < a[6];
    r[6] = temp + b[6];
    carry += r[6] < temp;

    temp = a[7] + carry;
    carry = temp < a[7];
    r[7] = temp + b[7];
    carry += r[7] < temp;

    return carry;
}

int sm2_bn_get_booth(const sm2_bn_type a, const unsigned int window_size, const unsigned int i) {
    const uint64_t mask = (1 << window_size) - 1;
    if (i == 0) {
        return (int) (a[0] << 1 & mask) - (int) (a[0] & mask);
    }

    unsigned int j = i * window_size - 1;
    const unsigned int n = j >> 6;
    j = j % 64;

    uint64_t w_bits = a[n] >> j;
    if (64 - j < (int) (window_size + 1) && n < 3) {
        w_bits |= a[n + 1] << (64 - j);
    }
    return (int) (w_bits & mask) - (int) ((w_bits >> 1) & mask);
}

void sm2_bn_from_hex(sm2_bn_type r, const char *hex) {
    uint8_t bytes[32];
    size_t len;
    hex2bytes(hex, 64, bytes, &len);
    sm2_bn_from_bytes(r, bytes);
}

int sm2_bn_printf(const char *label, const sm2_bn_type a) {
    sm2_printf("%s: %016llx%016llx%016llx%016llx\n", label, a[3], a[2], a[1], a[0]);
    return 1;
}


// GF(p)
// p = 2^256 - 2^224 - 2^96 + 2^64 - 1
//   = 0xfffffffeffffffffffffffffffffffffffffffff00000000ffffffffffffffff
const uint64_t SM2_BN_P[4] = {
    0xffffffffffffffff, 0xffffffff00000000, 0xffffffffffffffff, 0xfffffffeffffffff,
};

const uint64_t *sm2_bn_prime(void) {
    return &SM2_BN_P[0];
}

// 2^256 - p = 2^224 + 2^96 - 2^64 + 1
const uint64_t SM2_BN_NEG_P[4] = {
    1, ((uint64_t) 1 << 32) - 1, 0, ((uint64_t) 1 << 32),
};

void sm2_bn_mod_p_plus(sm2_bn_type r, const sm2_bn_type a, const sm2_bn_type b) {
    const uint64_t carry = sm2_bn_plus(r, a, b);
    if (carry) {
        // a + b - p = (a + b - 2^256) + (2^256 - p)
        sm2_bn_plus(r, r, SM2_BN_NEG_P);
        return;
    }
    if (sm2_bn_compare(r, SM2_BN_P) >= 0) {
        sm2_bn_subtraction(r, r, SM2_BN_P);
    }
}

void sm2_bn_mod_p_subtraction(sm2_bn_type r, const sm2_bn_type a, const sm2_bn_type b) {
    const uint64_t carry = sm2_bn_subtraction(r, a, b);
    if (carry) {
        // a - b + p = (a - b + 2^256) - (2^256 - p)
        sm2_bn_subtraction(r, r, SM2_BN_NEG_P);
    }
}

void sm2_bn_mod_p_double(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_mod_p_plus(r, a, a);
}

void sm2_bn_mod_p_triple(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_type t;
    sm2_bn_mod_p_plus(t, a, a);
    sm2_bn_mod_p_plus(r, t, a);
}

void sm2_bn_mod_p_negative(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_subtraction(r, SM2_BN_P, a);
}

void sm2_bn_mod_p_half(sm2_bn_type r, const sm2_bn_type a) {
    uint64_t carry = 0;
    if (a[0] & 1) {
        carry = sm2_bn_plus(r, a, SM2_BN_P);
    } else {
        sm2_bn_copy(r, a);
    }
    r[0] = (r[0] >> 1) | ((r[1] & 1) << 63);
    r[1] = (r[1] >> 1) | ((r[2] & 1) << 63);
    r[2] = (r[2] >> 1) | ((r[3] & 1) << 63);
    r[3] = (r[3] >> 1) | ((carry & 1) << 63);
}

// p' * p = -1 mod 2^256
// p' = -p^(-1) mod 2^256
//    = fffffffc00000001fffffffe00000000ffffffff000000010000000000000001
// sage: -(IntegerModRing(2^256)(p))^-1
const uint64_t SM2_BN_P_PRIME[4] = {
    0x0000000000000001, 0xffffffff00000001, 0xfffffffe00000000, 0xfffffffc00000001,
};


// mont(1) (mod p) = 2^256 mod p = 2^256 - p
const uint64_t *SM2_BN_MOD_P_MONT_ONE = SM2_BN_NEG_P;

// z = a*b
// c = (z + (z * p' mod 2^256) * p)/2^256
void sm2_bn_mod_p_mont_multiplication(sm2_bn_type r, const sm2_bn_type a, const sm2_bn_type b) {
    sm2_bbn_type z, t;
    // z = a * b
    sm2_bn_multiplication(z, a, b);

    // t = low(z) * p'
    sm2_bn_multiplication(t, z, SM2_BN_P_PRIME);

    // t = low(t) * p
    sm2_bn_multiplication(t, t, SM2_BN_P);

    // z = z + t
    const uint64_t carry = sm2_bbn_plus(z, z, t);

    // r = high(r)
    sm2_bn_copy(r, z + 4);

    if (carry) {
        sm2_bn_plus(r, r, SM2_BN_MOD_P_MONT_ONE);
    } else if (sm2_bn_compare(r, SM2_BN_P) >= 0) {
        sm2_bn_subtraction(r, r, SM2_BN_P);
    }
}

void sm2_bn_mod_p_mont_sqr(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_mod_p_mont_multiplication(r, a, a);
}

// mont(mont(a), 1) = aR * 1 * R^-1 (mod p) = a (mod p)
void sm2_bn_mod_p_from_mont(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_mod_p_mont_multiplication(r, a, SM2_BN_ONE);
}

// 2^512 (mod p)
const uint64_t SM2_BN_2e512mod_p[4] = {
    0x0000000200000003, 0x00000002ffffffff, 0x0000000100000001, 0x0000000400000002
};

// mont(a) = a * 2^256 (mod p) = mont_mul(a, 2^512 mod p)
void sm2_bn_mod_p_to_mont(const sm2_bn_type a, uint64_t r[4]) {
    sm2_bn_mod_p_mont_multiplication(r, a, SM2_BN_2e512mod_p);
}

// caller should check a != 0
void sm2_bn_mod_p_mont_inv(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_type a1, a2, a3, a4, a5;
    sm2_bn_mod_p_mont_sqr(a1, a);
    sm2_bn_mod_p_mont_multiplication(a2, a1, a);
    sm2_bn_mod_p_mont_sqr(a3, a2);
    sm2_bn_mod_p_mont_sqr(a3, a3);
    sm2_bn_mod_p_mont_multiplication(a3, a3, a2);
    sm2_bn_mod_p_mont_sqr(a4, a3);
    sm2_bn_mod_p_mont_sqr(a4, a4);
    sm2_bn_mod_p_mont_sqr(a4, a4);
    sm2_bn_mod_p_mont_sqr(a4, a4);
    sm2_bn_mod_p_mont_multiplication(a4, a4, a3);
    sm2_bn_mod_p_mont_sqr(a5, a4);
    for (int i = 1; i < 8; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a5, a5, a4);
    for (int i = 0; i < 8; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a5, a5, a4);
    for (int i = 0; i < 4; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a5, a5, a3);
    sm2_bn_mod_p_mont_sqr(a5, a5);
    sm2_bn_mod_p_mont_sqr(a5, a5);
    sm2_bn_mod_p_mont_multiplication(a5, a5, a2);
    sm2_bn_mod_p_mont_sqr(a5, a5);
    sm2_bn_mod_p_mont_multiplication(a5, a5, a);
    sm2_bn_mod_p_mont_sqr(a4, a5);
    sm2_bn_mod_p_mont_multiplication(a3, a4, a1);
    sm2_bn_mod_p_mont_sqr(a5, a4);
    for (int i = 1; i < 31; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a4, a5, a4);
    sm2_bn_mod_p_mont_sqr(a4, a4);
    sm2_bn_mod_p_mont_multiplication(a4, a4, a);
    sm2_bn_mod_p_mont_multiplication(a3, a4, a2);
    for (int i = 0; i < 33; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a2, a5, a3);
    sm2_bn_mod_p_mont_multiplication(a3, a2, a3);
    for (int i = 0; i < 32; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a2, a5, a3);
    sm2_bn_mod_p_mont_multiplication(a3, a2, a3);
    sm2_bn_mod_p_mont_multiplication(a4, a2, a4);
    for (int i = 0; i < 32; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a2, a5, a3);
    sm2_bn_mod_p_mont_multiplication(a3, a2, a3);
    sm2_bn_mod_p_mont_multiplication(a4, a2, a4);
    for (int i = 0; i < 32; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a2, a5, a3);
    sm2_bn_mod_p_mont_multiplication(a3, a2, a3);
    sm2_bn_mod_p_mont_multiplication(a4, a2, a4);
    for (int i = 0; i < 32; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(a2, a5, a3);
    sm2_bn_mod_p_mont_multiplication(a3, a2, a3);
    sm2_bn_mod_p_mont_multiplication(a4, a2, a4);
    for (int i = 0; i < 32; i++) {
        sm2_bn_mod_p_mont_sqr(a5, a5);
    }
    sm2_bn_mod_p_mont_multiplication(r, a4, a5);
}

// (p+1)/4 = 3fffffffbfffffffffffffffffffffffffffffffc00000004000000000000000
const uint64_t SM2_BN_SQRT_EXP[4] = {
    0x4000000000000000, 0xffffffffc0000000, 0xffffffffffffffff, 0x3fffffffbfffffff,
};

// n = 0xfffffffeffffffffffffffffffffffff7203df6b21c6052b53bbf40939d54123
const uint64_t SM2_BN_N[4] = {
    0x53bbf40939d54123, 0x7203df6b21c6052b, 0xffffffffffffffff, 0xfffffffeffffffff,
};

const uint64_t SM2_BN_N_MINUS_ONE[4] = {
    0x53bbf40939d54122, 0x7203df6b21c6052b, 0xffffffffffffffff, 0xfffffffeffffffff,
};

// 2^256 - n = 0x10000000000000000000000008dfc2094de39fad4ac440bf6c62abedd
const uint64_t SM2_BN_NEG_N[4] = {
    0xac440bf6c62abedd, 0x8dfc2094de39fad4, 0x0000000000000000, 0x0000000100000000,
};

void sm2_bn_mod_n_plus(sm2_bn_type r, const sm2_bn_type a, const sm2_bn_type b) {
    const uint64_t carry = sm2_bn_plus(r, a, b);

    if (carry) {
        // a + b - n = (a + b - 2^256) + (2^256 - n)
        sm2_bn_plus(r, r, SM2_BN_NEG_N);
        return;
    }

    if (sm2_bn_compare(r, SM2_BN_N) >= 0) {
        sm2_bn_subtraction(r, r, SM2_BN_N);
    }
}

void sm2_bn_mod_n_subtraction(sm2_bn_type r, const sm2_bn_type a, const sm2_bn_type b) {
    const uint64_t carry = sm2_bn_subtraction(r, a, b);
    if (carry) {
        // a - b + n = (a - b + 2^256) - (2^256 - n)
        sm2_bn_subtraction(r, r, SM2_BN_NEG_N);
    }
}


// n' = -n^(-1) mod 2^256
//    = 0x6f39132f82e4c7bc2b0068d3b08941d4df1e8d34fc8319a5327f9e8872350975
// sage: -(IntegerModRing(2^256)(n))^-1
const uint64_t SM2_BN_N_PRIME[4] = {
    0x327f9e8872350975, 0xdf1e8d34fc8319a5, 0x2b0068d3b08941d4, 0x6f39132f82e4c7bc,
};

const uint64_t *sm2_bn_order(void) {
    return &SM2_BN_N[0];
}

const uint64_t *sm2_bn_order_minus_one(void) {
    return &SM2_BN_N_MINUS_ONE[0];
}

// mont(1) (mod n) = 2^256 - n
const uint64_t *SM2_BN_MOD_N_MONT_ONE = SM2_BN_NEG_N;

void sm2_bn_mod_n_mont_multiplication(sm2_bn_type r, const sm2_bn_type a, const sm2_bn_type b) {
    sm2_bbn_type z, t;
    // z = a * b
    sm2_bn_multiplication(z, a, b);

    // t = low(z) * n'
    sm2_bn_multiplication(t, z, SM2_BN_N_PRIME);

    // t = low(t) * n
    sm2_bn_multiplication(t, t, SM2_BN_N);

    // z = z + t
    const uint64_t carry = sm2_bbn_plus(z, z, t);

    // r = high(r)
    sm2_bn_copy(r, z + 4);

    if (carry) {
        sm2_bn_plus(r, r, SM2_BN_MOD_N_MONT_ONE);
    } else if (sm2_bn_compare(r, SM2_BN_N) >= 0) {
        (void) sm2_bn_subtraction(r, r, SM2_BN_N);
    }
}

void sm2_bn_mod_n_mont_sqr(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_mod_n_mont_multiplication(r, a, a);
}

// n - 2 = 0xfffffffeffffffffffffffffffffffff7203df6b21c6052b53bbf40939d54121
const uint64_t SM2_BN_N_MINUS_TWO[4] = {
    0x53bbf40939d54121, 0x7203df6b21c6052b, 0xffffffffffffffff, 0xfffffffeffffffff,
};

void sm2_bn_mod_n_mont_inv(sm2_bn_type r, const sm2_bn_type a) {
    // expand sm2_bn_mod_n_mont_exp(r, a, SM2_BN_N_MINUS_TWO)
    sm2_bn_type t;
    sm2_bn_copy(t, a);

    for (int i = 0; i < 30; i++) {
        sm2_bn_mod_n_mont_sqr(t, t);
        sm2_bn_mod_n_mont_multiplication(t, t, a);
    }
    sm2_bn_mod_n_mont_sqr(t, t);
    for (int i = 0; i < 96; i++) {
        sm2_bn_mod_n_mont_sqr(t, t);
        sm2_bn_mod_n_mont_multiplication(t, t, a);
    }
    uint64_t w = SM2_BN_N_MINUS_TWO[1];
    for (int i = 0; i < 64; i++) {
        sm2_bn_mod_n_mont_sqr(t, t);
        if (w & 0x8000000000000000) {
            sm2_bn_mod_n_mont_multiplication(t, t, a);
        }
        w <<= 1;
    }
    w = SM2_BN_N_MINUS_TWO[0];
    for (int i = 0; i < 64; i++) {
        sm2_bn_mod_n_mont_sqr(t, t);
        if (w & 0x8000000000000000) {
            sm2_bn_mod_n_mont_multiplication(t, t, a);
        }
        w <<= 1;
    }
    sm2_bn_copy(r, t);
}

void sm2_bn_mod_n_inv(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_type mont_a;
    sm2_bn_mod_n_to_mont(a, mont_a);
    sm2_bn_mod_n_mont_inv(r, mont_a);
    sm2_bn_mod_n_from_mont(r, r);
}


// mont(mont(a), 1) = aR * 1 * R^-1 (mod n) = a (mod p)
void sm2_bn_mod_n_from_mont(sm2_bn_type r, const sm2_bn_type a) {
    sm2_bn_mod_n_mont_multiplication(r, a, SM2_BN_ONE);
}

// 2^512 (mod n) = 0x1eb5e412a22b3d3b620fc84c3affe0d43464504ade6fa2fa901192af7c114f20
const uint64_t SM2_BN_2e512mod_n[4] = {
    0x901192af7c114f20, 0x3464504ade6fa2fa, 0x620fc84c3affe0d4, 0x1eb5e412a22b3d3b,
};

// mont(a) = a * 2^256 (mod n) = mont_mul(a, 2^512 mod n)
void sm2_bn_mod_n_to_mont(const sm2_bn_type a, uint64_t r[4]) {
    sm2_bn_mod_n_mont_multiplication(r, a, SM2_BN_2e512mod_n);
}


// Jacobian Point with Montgomery coordinates
void sm2_bn_point_set_infinity(sm2_bn_point_type *P) {
    sm2_bn_copy(P->X, SM2_BN_MOD_P_MONT_ONE);
    sm2_bn_copy(P->Y, SM2_BN_MOD_P_MONT_ONE);
    sm2_bn_set_zero(P->Z);
}

// point at infinity should be like (k^2 : k^3 : 0), k in [0, p-1]
int sm2_bn_point_is_at_infinity(const sm2_bn_point_type *P) {
    if (sm2_bn_is_zero(P->Z)) {
        sm2_bn_type X_cub, Y_sqr;
        sm2_bn_mod_p_mont_sqr(X_cub, P->X);
        sm2_bn_mod_p_mont_multiplication(X_cub, X_cub, P->X);
        sm2_bn_mod_p_mont_sqr(Y_sqr, P->Y);
        if (sm2_bn_compare(X_cub, Y_sqr) != 0) {
            print_log();
            return 0;
        }
        return 1;
    }
    return 0;
}

// mont(b), b = 0x28e9fa9e9d9f5e344d5a9e4bcf6509a7f39789f515ab8f92ddbcbd414d940e93
const uint64_t SM2_BN_MOD_P_MONT_B[4] = {
    0x90d230632bc0dd42, 0x71cf379ae9b537ab, 0x527981505ea51c3c, 0x240fe188ba20e2c8,
};

int sm2_bn_point_is_on_curve(const sm2_bn_point_type *P) {
    sm2_bn_type t0, t1;
    if (sm2_bn_compare(P->Z, SM2_BN_MOD_P_MONT_ONE) == 0) {
        // if Z == 1, check y^2 + 3*x == x^3 + b
        sm2_bn_mod_p_mont_sqr(t0, P->Y);
        sm2_bn_mod_p_plus(t0, t0, P->X);
        sm2_bn_mod_p_plus(t0, t0, P->X);
        sm2_bn_mod_p_plus(t0, t0, P->X);
        sm2_bn_mod_p_mont_sqr(t1, P->X);
        sm2_bn_mod_p_mont_multiplication(t1, t1, P->X);
        sm2_bn_mod_p_plus(t1, t1, SM2_BN_MOD_P_MONT_B);
    } else {
        sm2_bn_type t2;
        // check Y^2 + 3 * X * Z^4 == X^3 + b * Z^6
        // if Z == 0, Y^2 == X^3, i.e. Y == X is checked
        sm2_bn_mod_p_mont_sqr(t0, P->Y);
        sm2_bn_mod_p_mont_sqr(t1, P->Z);
        sm2_bn_mod_p_mont_sqr(t2, t1);
        sm2_bn_mod_p_mont_multiplication(t1, t1, t2);
        sm2_bn_mod_p_mont_multiplication(t1, t1, SM2_BN_MOD_P_MONT_B);
        sm2_bn_mod_p_mont_multiplication(t2, t2, P->X);
        sm2_bn_mod_p_plus(t0, t0, t2);
        sm2_bn_mod_p_plus(t0, t0, t2);
        sm2_bn_mod_p_plus(t0, t0, t2);
        sm2_bn_mod_p_mont_sqr(t2, P->X);
        sm2_bn_mod_p_mont_multiplication(t2, t2, P->X);
        sm2_bn_mod_p_plus(t1, t1, t2);
    }
    if (sm2_bn_compare(t0, t1) != 0) {
        return 0;
    }
    return 1;
}

int sm2_bn_point_get_xy(const sm2_bn_point_type *P, uint64_t x[4], uint64_t y[4]) {
    if (sm2_bn_point_is_at_infinity(P) == 1) {
        sm2_bn_set_zero(x);
        if (y) {
            sm2_bn_set_zero(y);
        }
        return 0;
    }
    if (sm2_bn_compare(P->Z, SM2_BN_MOD_P_MONT_ONE) == 0) {
        sm2_bn_mod_p_from_mont(x, P->X);
        if (y) {
            sm2_bn_mod_p_from_mont(y, P->Y);
        }
    } else {
        sm2_bn_type z_inv;
        sm2_bn_mod_p_mont_inv(z_inv, P->Z);
        if (y) {
            sm2_bn_mod_p_mont_multiplication(y, P->Y, z_inv);
        }
        sm2_bn_mod_p_mont_sqr(z_inv, z_inv);
        sm2_bn_mod_p_mont_multiplication(x, P->X, z_inv);
        sm2_bn_mod_p_from_mont(x, x);
        if (y) {
            sm2_bn_mod_p_mont_multiplication(y, y, z_inv);
            sm2_bn_mod_p_from_mont(y, y);
        }
    }
    return 1;
}

void sm2_bn_point_double(sm2_bn_point_type *R, const sm2_bn_point_type *A) {
    const uint64_t *X1 = A->X, *Y1 = A->Y, *Z1 = A->Z;
    uint64_t *X3 = R->X, *Y3 = R->Y, *Z3 = R->Z;
    sm2_bn_type S, M, Z_sqr, tmp;

    // 1. S = 2Y
    sm2_bn_mod_p_double(S, Y1);

    // 2. Z_sqr = Z^2
    sm2_bn_mod_p_mont_sqr(Z_sqr, Z1);

    // 3. S = S^2 = 4Y^2
    sm2_bn_mod_p_mont_sqr(S, S);

    // 4. Z = Z*Y
    sm2_bn_mod_p_mont_multiplication(Z3, Z1, Y1);

    // 5. Z = 2*Z = 2*Y*Z ===> Z3
    sm2_bn_mod_p_double(Z3, Z3);

    // 6. M = X + Z_sqr = X + Z^2
    sm2_bn_mod_p_plus(M, X1, Z_sqr);

    // 7. Z_sqr = X - Z_sqr = X - Z^2
    sm2_bn_mod_p_subtraction(Z_sqr, X1, Z_sqr);

    // 8. Y = S^2 = 16Y^4
    sm2_bn_mod_p_mont_sqr(Y3, S);

    // 9. Y = Y/2 = 8Y^4
    sm2_bn_mod_p_half(Y3, Y3);

    // 10. M = M * Z_sqr = (X + Z^2)*(X - Z^2) = X^2 - Z^4
    sm2_bn_mod_p_mont_multiplication(M, M, Z_sqr);

    // 11. M = 3M = 3X^2 - 3Z^4
    sm2_bn_mod_p_triple(M, M);

    // 12. S = S * X = 4X*Y^2
    sm2_bn_mod_p_mont_multiplication(S, S, X1);

    // 13. tmp = 2 * S = 8X*Y^2
    sm2_bn_mod_p_double(tmp, S);

    // 14. X = M^2 = (3X^2 - 3Z^4)^2
    sm2_bn_mod_p_mont_sqr(X3, M);

    // 15. X = X - tmp = (3X^2 - 3Z^4)^2 - 8X*Y^2 ===> X3
    sm2_bn_mod_p_subtraction(X3, X3, tmp);

    // 16. S = S - X3 = 4X*Y^2 - X3
    sm2_bn_mod_p_subtraction(S, S, X3);

    // 17. S = S * M = (3X^2 - 3Z^4)*(4X*Y^2 - X3)
    sm2_bn_mod_p_mont_multiplication(S, S, M);

    // 18. Y = S - Y = (3X^2 - 3Z^4)*(4X*Y^2 - X3) - 8Y^4 ===> Y3
    sm2_bn_mod_p_subtraction(Y3, S, Y3);
}

/*
  (X1:Y1:Z1) + (X2:Y2:Z2) => (X3:Y3:Z3)

	A = Y2 * Z1^3 - Y1 * Z2^3
	B = X2 * Z1^2 - X1 * Z2^2

	X3 = A^2 - B^2 * (X2 * Z1^2 + X1 * Z2^2)
	   = A^2 - B^3 - 2 * B^2 * X1 * Z2^2
	Y3 = A * (X1 * B^2 * Z2^2 - X3) - Y1 * B^3 * Z2^3
	Z3 = B * Z1 * Z2

  P + (-P) = (X:Y:Z) + (k^2*X : k^3*Y : k*Z) => (0:0:0)
*/
void sm2_bn_point_plus(sm2_bn_point_type *r, const sm2_bn_point_type *a, const sm2_bn_point_type *b) {
    sm2_bn_type U1, S1, U2, S2;
    sm2_bn_type Z1_sqr, Z2_sqr;
    sm2_bn_type H, R;
    sm2_bn_type H_sqr, R_sqr, H_cub;

    sm2_bn_type res_x, res_y, res_z;

    const uint64_t *in1_x = a->X, *in1_y = a->Y, *in1_z = a->Z;
    const uint64_t *in2_x = b->X, *in2_y = b->Y, *in2_z = b->Z;

    /*
    * Infinity in encoded as (x,y,0)
    */
    uint64_t in1_infinity = (in1_z[0] | in1_z[1] | in1_z[2] | in1_z[3]);
    uint64_t in2_infinity = (in2_z[0] | in2_z[1] | in2_z[2] | in2_z[3]);

    in1_infinity = uint64_t_is_zero(in1_infinity);
    in2_infinity = uint64_t_is_zero(in2_infinity);

    sm2_bn_mod_p_mont_sqr(Z2_sqr, in2_z); /* Z2^2 */
    sm2_bn_mod_p_mont_sqr(Z1_sqr, in1_z); /* Z1^2 */

    sm2_bn_mod_p_mont_multiplication(S1, Z2_sqr, in2_z); /* S1 = Z2^3 */
    sm2_bn_mod_p_mont_multiplication(S2, Z1_sqr, in1_z); /* S2 = Z1^3 */

    sm2_bn_mod_p_mont_multiplication(S1, S1, in1_y); /* S1 = Y1*Z2^3 */
    sm2_bn_mod_p_mont_multiplication(S2, S2, in2_y); /* S2 = Y2*Z1^3 */
    sm2_bn_mod_p_subtraction(R, S2, S1); /* R = S2 - S1 */

    sm2_bn_mod_p_mont_multiplication(U1, in1_x, Z2_sqr); /* U1 = X1*Z2^2 */
    sm2_bn_mod_p_mont_multiplication(U2, in2_x, Z1_sqr); /* U2 = X2*Z1^2 */
    sm2_bn_mod_p_subtraction(H, U2, U1); /* H = U2 - U1 */

    /*
    * This should not happen during sign/ecdh, so no constant time violation
    */
    if (sm2_bn_is_equal(U1, U2) && !in1_infinity && !in2_infinity) {
        if (sm2_bn_is_equal(S1, S2)) {
            sm2_bn_point_double(r, a);
            return;
        }
        memset(r, 0, sizeof(*r));
        return;
    }

    sm2_bn_mod_p_mont_sqr(R_sqr, R); /* R^2 */
    sm2_bn_mod_p_mont_multiplication(res_z, H, in1_z); /* Z3 = H*Z1*Z2 */

    sm2_bn_mod_p_mont_sqr(H_sqr, H); /* H^2 */
    sm2_bn_mod_p_mont_multiplication(res_z, res_z, in2_z); /* Z3 = H*Z1*Z2 */

    sm2_bn_mod_p_mont_multiplication(H_cub, H_sqr, H); /* H^3 */
    sm2_bn_mod_p_mont_multiplication(U2, U1, H_sqr); /* U1*H^2 */

    sm2_bn_mod_p_double(H_sqr, U2); /* 2*U1*H^2 */

    sm2_bn_mod_p_subtraction(res_x, R_sqr, H_sqr);
    sm2_bn_mod_p_subtraction(res_x, res_x, H_cub);

    sm2_bn_mod_p_subtraction(res_y, U2, res_x);

    sm2_bn_mod_p_mont_multiplication(S2, S1, H_cub);
    sm2_bn_mod_p_mont_multiplication(res_y, R, res_y);

    sm2_bn_mod_p_subtraction(res_y, res_y, S2);

    sm2_bn_copy_conditional(res_x, in2_x, in1_infinity);
    sm2_bn_copy_conditional(res_y, in2_y, in1_infinity);
    sm2_bn_copy_conditional(res_z, in2_z, in1_infinity);

    sm2_bn_copy_conditional(res_x, in1_x, in2_infinity);
    sm2_bn_copy_conditional(res_y, in1_y, in2_infinity);
    sm2_bn_copy_conditional(res_z, in1_z, in2_infinity);

    memcpy(r->X, res_x, sizeof(res_x));
    memcpy(r->Y, res_y, sizeof(res_y));
    memcpy(r->Z, res_z, sizeof(res_z));
}

void sm2_bn_point_negative(sm2_bn_point_type *R, const sm2_bn_point_type *P) {
    sm2_bn_copy(R->X, P->X);
    sm2_bn_mod_p_negative(R->Y, P->Y);
    sm2_bn_copy(R->Z, P->Z);
}

void sm2_bn_point_subtraction(sm2_bn_point_type *R, const sm2_bn_point_type *A, const sm2_bn_point_type *B) {
    sm2_bn_point_type neg_B;
    sm2_bn_point_negative(&neg_B, B);
    sm2_bn_point_plus(R, A, &neg_B);
}

void sm2_bn_point_mul_pre_compute(const sm2_bn_point_type *P, sm2_bn_point_type T[16]) {
    memcpy(&T[0], P, sizeof(sm2_bn_point_type));
    if (sm2_bn_is_equal(P->Z, SM2_BN_MOD_P_MONT_ONE) == 1) {
        const SM2_BN_AFFINE_POINT *P_ = (const SM2_BN_AFFINE_POINT *) P;
        sm2_bn_point_double(&T[1], &T[0]);
        sm2_bn_point_add_affine(&T[2], &T[1], P_);
        sm2_bn_point_double(&T[3], &T[1]);
        sm2_bn_point_add_affine(&T[4], &T[3], P_);
        sm2_bn_point_double(&T[5], &T[2]);
        sm2_bn_point_add_affine(&T[6], &T[5], P_);
        sm2_bn_point_double(&T[7], &T[3]);
        sm2_bn_point_add_affine(&T[8], &T[7], P_);
        sm2_bn_point_double(&T[9], &T[4]);
        sm2_bn_point_add_affine(&T[10], &T[9], P_);
        sm2_bn_point_double(&T[11], &T[5]);
        sm2_bn_point_add_affine(&T[12], &T[11], P_);
        sm2_bn_point_double(&T[13], &T[6]);
        sm2_bn_point_add_affine(&T[14], &T[13], P_);
        sm2_bn_point_double(&T[15], &T[7]);
    } else {
        sm2_bn_point_double(&T[2 - 1], &T[0]);
        sm2_bn_point_double(&T[4 - 1], &T[2 - 1]);
        sm2_bn_point_double(&T[8 - 1], &T[4 - 1]);
        sm2_bn_point_double(&T[16 - 1], &T[8 - 1]);
        sm2_bn_point_plus(&T[3 - 1], &T[2 - 1], P);
        sm2_bn_point_double(&T[6 - 1], &T[3 - 1]);
        sm2_bn_point_double(&T[12 - 1], &T[6 - 1]);
        sm2_bn_point_plus(&T[5 - 1], &T[3 - 1], &T[2 - 1]);
        sm2_bn_point_double(&T[10 - 1], &T[5 - 1]);
        sm2_bn_point_plus(&T[7 - 1], &T[4 - 1], &T[3 - 1]);
        sm2_bn_point_double(&T[14 - 1], &T[7 - 1]);
        sm2_bn_point_plus(&T[9 - 1], &T[4 - 1], &T[5 - 1]);
        sm2_bn_point_plus(&T[11 - 1], &T[6 - 1], &T[5 - 1]);
        sm2_bn_point_plus(&T[13 - 1], &T[7 - 1], &T[6 - 1]);
        sm2_bn_point_plus(&T[15 - 1], &T[8 - 1], &T[7 - 1]);
    }
}

void sm2_bn_point_mul_ex(sm2_bn_point_type *R, const uint64_t k[4], const sm2_bn_point_type *P_table) {
    const int window_size = 5;
    const int n = (256 + window_size - 1) / window_size;
    int R_infinity = 1;
    for (int i = n - 1; i >= 0; i--) {
        const int booth = sm2_bn_get_booth(k, window_size, i);
        if (R_infinity) {
            if (booth != 0) {
                *R = P_table[booth - 1];
                R_infinity = 0;
            }
        } else {
            sm2_bn_point_double(R, R);
            sm2_bn_point_double(R, R);
            sm2_bn_point_double(R, R);
            sm2_bn_point_double(R, R);
            sm2_bn_point_double(R, R);

            if (booth > 0) {
                sm2_bn_point_plus(R, R, &P_table[booth - 1]);
            } else if (booth < 0) {
                sm2_bn_point_subtraction(R, R, &P_table[-booth - 1]);
            }
        }
    }
    if (R_infinity) {
        memset(R, 0, sizeof(*R));
    }
}


int sm2_bn_point_printf(const char *label, const sm2_bn_point_type *P) {
    if (sm2_bn_point_is_at_infinity(P) == 1) {
        sm2_printf("%s: point_at_infinity\n", label);
    } else {
        uint8_t bytes[64];
        sm2_bn_point_to_bytes(P, bytes);
        sm2_printf_bytes(label, bytes, 64);
    }
    return 1;
}

void sm2_bn_point_copy_affine(sm2_bn_point_type *R, const SM2_BN_AFFINE_POINT *P) {
    memcpy(R, P, sizeof(SM2_BN_AFFINE_POINT));
    sm2_bn_copy(R->Z, SM2_BN_MOD_P_MONT_ONE);
}

void sm2_bn_point_add_affine(sm2_bn_point_type *r, const sm2_bn_point_type *a, const SM2_BN_AFFINE_POINT *b) {
    sm2_bn_type U2, S2, Z1sqr;
    sm2_bn_type H, R;
    sm2_bn_type H_sqr, R_sqr, H_cub;

    sm2_bn_type res_x, res_y, res_z;

    const uint64_t *in1_x = a->X, *in1_y = a->Y, *in1_z = a->Z;
    const uint64_t *in2_x = b->x, *in2_y = b->y;

    /*
    * Infinity in encoded as (x,y,0)
    */
    uint64_t in1_infinity = (in1_z[0] | in1_z[1] | in1_z[2] | in1_z[3]);

    /*
    * In affine representation we encode infinity as (0,0), which is
    * not on the curve, so it is OK
    */
    uint64_t in2_infinity = (in2_x[0] | in2_x[1] | in2_x[2] | in2_x[3] | in2_y[0] | in2_y[1] | in2_y[2] | in2_y[3]);

    in1_infinity = uint64_t_is_zero(in1_infinity);
    in2_infinity = uint64_t_is_zero(in2_infinity);


    /* Z1^2 */
    sm2_bn_mod_p_mont_sqr(Z1sqr, in1_z);

    /* U2 = X2*Z1^2 */
    sm2_bn_mod_p_mont_multiplication(U2, in2_x, Z1sqr);
    /* H = U2 - U1 */
    sm2_bn_mod_p_subtraction(H, U2, in1_x);
    sm2_bn_mod_p_mont_multiplication(S2, Z1sqr, in1_z); /* S2 = Z1^3 */
    sm2_bn_mod_p_mont_multiplication(res_z, H, in1_z); /* Z3 = H*Z1*Z2 */
    sm2_bn_mod_p_mont_multiplication(S2, S2, in2_y); /* S2 = Y2*Z1^3 */
    sm2_bn_mod_p_subtraction(R, S2, in1_y); /* R = S2 - S1 */
    sm2_bn_mod_p_mont_sqr(H_sqr, H); /* H^2 */
    sm2_bn_mod_p_mont_sqr(R_sqr, R); /* R^2 */
    sm2_bn_mod_p_mont_multiplication(H_cub, H_sqr, H); /* H^3 */
    sm2_bn_mod_p_mont_multiplication(U2, in1_x, H_sqr); /* U1*H^2 */
    sm2_bn_mod_p_double(H_sqr, U2); /* 2*U1*H^2 */
    sm2_bn_mod_p_subtraction(res_x, R_sqr, H_sqr);
    sm2_bn_mod_p_subtraction(res_x, res_x, H_cub);
    sm2_bn_mod_p_subtraction(H, U2, res_x);
    sm2_bn_mod_p_mont_multiplication(S2, in1_y, H_cub);
    sm2_bn_mod_p_mont_multiplication(H, H, R);
    sm2_bn_mod_p_subtraction(res_y, H, S2);
    sm2_bn_copy_conditional(res_x, in2_x, in1_infinity);
    sm2_bn_copy_conditional(res_x, in1_x, in2_infinity);
    sm2_bn_copy_conditional(res_y, in2_y, in1_infinity);
    sm2_bn_copy_conditional(res_y, in1_y, in2_infinity);
    sm2_bn_copy_conditional(res_z, SM2_BN_MOD_P_MONT_ONE, in1_infinity);
    sm2_bn_copy_conditional(res_z, in1_z, in2_infinity);

    memcpy(r->X, res_x, sizeof(res_x));
    memcpy(r->Y, res_y, sizeof(res_y));
    memcpy(r->Z, res_z, sizeof(res_z));
}

void sm2_bn_point_sub_affine(sm2_bn_point_type *R, const sm2_bn_point_type *A, const SM2_BN_AFFINE_POINT *B) {
    SM2_BN_AFFINE_POINT neg_B;
    sm2_bn_copy(neg_B.x, B->x);
    sm2_bn_mod_p_negative(neg_B.y, B->y);
    sm2_bn_point_add_affine(R, A, &neg_B);
}

extern const uint64_t sm2_bn_pre_compute_table[37][64 * 4 * 2];
static SM2_BN_AFFINE_POINT (*g_pre_comp)[64] = (SM2_BN_AFFINE_POINT (*)[64]) sm2_bn_pre_compute_table;

void sm2_bn_point_mul_generator(sm2_bn_point_type *R, const sm2_bn_type k) {
    const unsigned int window_size = 7;
    int R_infinity = 1;
    const unsigned int n = (256 + window_size - 1) / window_size;

    for (int i = (int) n - 1; i >= 0; i--) {
        const int booth = sm2_bn_get_booth(k, window_size, i);
        if (R_infinity) {
            if (booth != 0) {
                sm2_bn_point_copy_affine(R, &g_pre_comp[i][booth - 1]);
                R_infinity = 0;
            }
        } else {
            if (booth > 0) {
                sm2_bn_point_add_affine(R, R, &g_pre_comp[i][booth - 1]);
            } else if (booth < 0) {
                sm2_bn_point_sub_affine(R, R, &g_pre_comp[i][-booth - 1]);
            }
        }
    }
    if (R_infinity) {
        sm2_bn_point_set_infinity(R);
    }
}

// point_at_infinity can not be encoded/decoded to/from bytes
int sm2_bn_point_from_bytes(sm2_bn_point_type *P, const uint8_t in[64]) {
    sm2_bn_from_bytes(P->X, in);
    if (sm2_bn_compare(P->X, sm2_bn_prime()) >= 0) {
        print_log();
        return -1;
    }
    sm2_bn_from_bytes(P->Y, in + 32);
    if (sm2_bn_compare(P->Y, sm2_bn_prime()) >= 0) {
        print_log();
        return -1;
    }

    // point_at_infinity
    if (sm2_bn_is_zero(P->X) == 1 && sm2_bn_is_zero(P->Y) == 1) {
        sm2_bn_point_set_infinity(P);
        return 0;
    }

    sm2_bn_mod_p_to_mont(P->X, P->X);
    sm2_bn_mod_p_to_mont(P->Y, P->Y);
    sm2_bn_copy(P->Z, SM2_BN_MOD_P_MONT_ONE);

    if (sm2_bn_point_is_on_curve(P) != 1) {
        print_log();
        return -1;
    }
    return 1;
}

int sm2_bn_point_from_hex(sm2_bn_point_type *P, const char *hex) {
    uint8_t bytes[64];
    size_t len;
    int ret;
    hex2bytes(hex, 128, bytes, &len);
    if ((ret = sm2_bn_point_from_bytes(P, bytes)) < 0) {
        print_log();
        return -1;
    }
    return ret;
}

// point_at_infinity should not to_bytes
int sm2_bn_point_to_bytes(const sm2_bn_point_type *P, uint8_t out[64]) {
    sm2_bn_type x, y;
    const int ret = sm2_bn_point_get_xy(P, x, y);
    sm2_bn_to_bytes(x, out);
    sm2_bn_to_bytes(y, out + 32);
    return ret;
}
